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We present and implement an efficient variational method to simulate two-dimensional finite size 
fermionic quantum systems by fermionic projected entangled pair states. The approach differs from 
the original one due to the fact that there is no need for an extra string-bond for contracting the 
tensor network. The method is tested on a bi-linear fermionic model on a square lattice for sizes up 
to ten by ten where good relative accuracy is achieved. Qualitatively good results are also obtained 
for an interacting fermionic system. 
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I. INTRODUCTION 

The theoretical study of quantum many-body systems 
presents one of the most challenging tasks of condensed 
matter physics, computational physics and quantum 
chemistry. Several approaches have been proposed to 
study quantum many-body systems e.g. quantum monte 
carlo (QMC), dynamical mean field theory (DMFT) , 
density-matrix renormalization group (DMRG) 1 2 /tensor 
network methodflH. The latter are best suited to de- 
scribe physical systems of local Hamiltonians at zero 
temperature in one spatial dimension for which it was 
showrPHU that they can be well approximated by matrix 
product states (MPS). 

A generalization of MPS algorithms to two spatial di- 
mensions was given by projected entangled pair states 
(PEPS)^ 9 algorithms where the quantum state is de- 
scribed in terms of entangled pairs on lattice bonds. 
Those states capture the entanglement structure needed 
to represent states that obey an area law^Sl, and there are 
strong arguments why every ground state of a gapped 
two-dimensional local Hamiltonian can be efficiently rep- 
resented as a 

All these methods were originally constructed to sim- 
ulate quantum spin systems whereas practical problems 
in condensed matter physics and quantum chemistry are 
more often of fermionic nature. A notorious example 
is the Fermi-Hubbard model which is believed to be a 
good candidate for high-temperature superconductivity. 
For local Hamiltonians in one spatial dimension, the dis- 
tinction between spins and fermions is irrelevant as any 
physical fermionic model can be transformed by Jordan- 
Wigner transformation to a spin model where the locality 
of interaction is preserved. This is not the case in two- 
dimensional systems where such transformation would in 
general convert local interactions to non-local strings op- 
erators. In the case of the ladders it is in principle possi- 
ble to use linear DMRG methods if all the symmetries are 
exploited^ but such an approach is clearly not scalable. 

Around the same time, two independent approaches 
to simulate two-dimensional fermionic systems were pro- 
posed: a generalization of the multiscale entanglement 
renormalization (MERA)^ to fermionic systems^ and 
the description of two-dimensional fermionic systems in 



terms of fermionic projected pair states (fPEPS^l In 
the latter, PEPS were generalized to fermionic systems 
in a natural way by considering entangled fermionic 
pairs instead of entangled spin pairs. It was also shown 
that this is a good ansatz and is in principle able to 
parametrize ground states of gapped fermionic models. 
However, in the case of the fPEPS the sign problem was 
not solved in a completely satisfactory way as there was 
still a need for increasing the bond dimension with a fac- 
tor of two if one were to contract this fermionic PEPS us- 
ing the standard procedure. The efficient contraction of 
fermionic MERA was reinterpreted in terms of fermionic 
swap and jump rule d 15 ! 1 ^ which allowed to generalize the 
sign-free contraction also to arbitrary fermionic tensor 
networks, as first reported in RefP^ and subsequently in 
RefP. In these papers, it was also sketched how this for- 
malism can be used to contract general fPEPS. The first 
fPEPS simulations albeit without the s ign-free contrac- 
tion rules, were performed in Refs.LU*L!3D under the name 
Graded PEPS, and very promising numerical results were 
reported. Finally, the full sign-free fPEPS algorithm for 
infinite lattices was implemented, together with interest- 
ing numerical results on interacting fermions and the t-J 
model, in Refill where also an explicit scheme for con- 
tracting finite-size fPEPS was given. A crucial element 
in all those approaches was the realization of a simple 
fermionic swapping rule which will also play an impor- 
tant role in this paper where we focus on construction of 
the finite-size fermionic PEPS algorithm. 

The most obvious advantage of the finite lattice PEPS 
over the infinite PEPS (iPEPS) algorithm is that no as- 
sumption of translation invariance symmetry is required. 
The finite size PEPS algorithm is therefore well suited 
to simulate physical systems with an unknown transla- 
tion invariance pattern or translation non-invariant (dis- 
ordered or noisy) systems. The only input information 
for the finite PEPS algorithm is the Hamiltonian opera- 
tor and, possibly, the parity of the ground state. How- 
ever, if the symmetries are known, they can be embedded 
natural! jE3. 

In this paper, we address the fermionic PEPS 11 en- 
tirely in terms of fermions without introducing any addi- 
tional bonds between lattice sites but rather embedding 
all fermionic signs locally. The crucial element in such de- 
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scription is a fe rmionic rule used to swap two fermionic 
tensor operators^ 1121221 . This way the complexity of the 
method exactly translates to the conventional PEPS for 
quantum spin systems (strictly speaking, it is even more 
efficient due to the parity constraints). We are able to 
efficiently calculate expectation values of arbitrary oper- 
ators for a given fPEPS state and efficiently optimize the 
ground state approximation for an arbitrary fermionic 
system on a rectangular lattice. We test the method on 
an integrable quadratic model on a square lattice and 
compare the ground state energy and the total particle 
number to the exact values. 



II. FERMIONIC PROJECTED ENTANGLED 
PAIR STATES 

We start with rewriting the orig inal fPEPS ansatzSS m 
an alternative way which will allow fermionic manipula- 
tions and consequently sign-free contraction of fermionic 
PEPS states. As given in 11 , a quantum state of a 
fermionic system on a square lattice can be described 
in terms of fermionic entangled pair states as 

l*> = w„ tla ,s n Qt* II //,, II vy o) (i) 

i,j i,j i,3 

where to each site (i,j) four auxiliary fermions are as- 
sociated: on,ji Pi.ji li.j arL d 5i t j, connecting the site 
to the respective left, right, upper and lower neigh- 
bor. The entangled pairs on the horizontal and vertical 
bonds are created (up to the normalization) by operators 

H itj = and V itj = l+^'Ti+i.j: respectively, 

and the projection to the space of physical fermions is 
given by projectors 

Qi,j — ^lrudk C \,j a i,jfii,jli,j&i,j- (2) 

For brevity, we omit the summation symbol where it 
is understood that the summation takes place over 
all indices that appear both in sub-script and super- 
script. The expression is traced over the space of virtual 
fermions which is formally designated by the operator 
Wa,/3,7,<5 which mimics the vacuum of virtual particles 
in the subscript, e.g. = Y\ u a v a t an d Wa,p,-y,s = 

WgWpWjWs. Note that W & = We will use~a~dash 
notation when referring to sequences m = (mi, m2, ■ ■ ■) 
or tensors A m = A mi . m2j ... with the rank given by the 
context. 

A fundamental feature of fermionic systems due to 
the causality is that the system is always in a state 
with a well defined parity, (W) (*|c w |*) = 0. When 
the ansatz is used to describe the ground state of 
a fermionic system, one can therefore assume that the 
projection operators are either parity preserving 

{Pi.j = 0) or parity violating (Pi,j = 1) and can be de- 
scribed by only half of the tensor elements of A} 1 '-'' , i.e. 



(l+r + u+d+k) mod 2 ^ P id A^ dk = 0. The practi- 
cal consequence of such assumption is that the projection 
operators ^ either commute or anti-commutc. 

Let us consider a mxn lattice of fermions and choose 
the row-major order in ([I]) by multiplying the projection 
operators Q v by the entanglement creation operators H v 
and V v in ([I]). This results in a description in terms of 
two types of virtual fermions a v and 7„ on horizontal and 
vertical bonds, respectively, 

I*) = W &1 A m , n ■ ■ ■ A hn ■ ■ • Ai,i|0) (3) 

with in general non-commuting operators 

A ■ - A [i,j] r n a l n U -y u ^ d f41 

of the same parity as the corresponding Qij, i.e. ei- 
ther parity preserving (Pij — 0) or swapping (Pij = 1). 
Again, we use the operator W an to mimic the contrac- 
tion over virtual particles. 

A. Expectation values 

A starting point in the computation with fermionic 
tensor product states is the calculation of expectation 
values of arbitrary operators. Due to the linearity it is 
sufficient to calculate the expectation value of an arbi- 
trary product operator 

O = O m . n ■ ■ ■ O x , n ■ ■ • O1.1 (5) 

where Ojj are single-site operators of a well defined par- 
ity pi j. Explicitly, each can be either parity pre- 
serving {pi.j — 0) in which case it can be written as 

= Oofi°i,j c i j + ^1,1 c lj c ij f° r some coefficients 

Oqq' and Oil\ or P ar ity swapping (pij = 1) such as 

Oi,j — O^'fcij +Oi ) o' c lj f° r some coefficients Oq'i and 

0[ ■ The expectation value is formally written as 

(#|0|*) = (0\A'li---A'l, n OA m , n ---Ai,i\0) (6) 

where the conjugated state is described by a comple- 
mentary set of virtual fermions designated by a prime, 
i.e. 

Exact contraction of such a tensor network, albeit pos- 
sible, is inefficient due to the contraction order specified 
in ([6]) . In order to contract the fermionic tensor network 
efficiently, one must be able to first contract over the 
physical modes and then contract the double layer in an 
approximate way^. In both steps one must be able to 
swap the contraction order between two tensor operators 
sharing a common contraction leg which, as we will show, 
is possible due to the parity constraints in fermionic ten- 
sor network. The latter step is performed by merging 
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FIG. 1. Swapping of two operators AB — BA. 



rows together and representing the double row by a sin- 
gle row, such that the horizontal bond dimensions remain 
finite. In order to contract over the physical modes in the 
first step, the tensor network must be written in a way 
where both tensors corresponding to a specific site appear 
together, such as A\jOi t jAi t j which are of a globally de- 
fined parity pij regardless of parity Pjj of A(j. Indeed, 
substituting the result by effective operators A u and 
using a fermionic rule explained in the following, that is 
exactly what we are able to achieve. Let us in the fol- 
lowing present a rule which will allow us to rewrite e.g. 
A'\ jA'l j +1 — A'i j^A'i j which is needed to reverse the 
contraction order in the conj ugate la yer. An equivalent 
rule was already used in Refe r 5 * 1 ' 20 \ 

Fermionic swap rule: let us define an arbitrary opera- 
tor A and an operator B of well defined parity pb as 



a a 1 1 a T 

A = Aiara^-a-^ 



and 



B = B L t,rb{ L ^% (8) 



where we use a notation c— = c™ 1 c™ 2 ■ ■ • where Cj rep- 
resents fermionic annihilation operators. Note that any 
superposition of products of fermionic operators can be 
written in this form. Then the following statement can 
be made: t'he product AB contracted over all common 
modes, here explicitly denoted by 7 = (71,72,...), can 
be written in a reverse order ( ? ) as 

W*ABW = W^BAW with W = = ]J Ijl) (9) 



where A and B are obtained independently from A and 
B, respectively, in addition to a global parity sign ps, as 



.4 
B 



(10) 



where (-1)-= (-l) a; i+^+-. 

The fermionic swap rule can be proven in a straight- 
forward way by writing the above ansatz and reordering 
the fermionic operators where all fermionic signs cancel 
except for the parity sign of the contraction product op- 
erator. It is however crucial to assume that at least one 
of the operators is of well defined parity. Otherwise, such 
decomposition is impossible since (— l) xy cannot be de- 
composed to a product f{x)g{y) for {0, 1} and any 
functions / and g. 

In the context of this paper, both A and B will always 
be of well defined parity. In such a case, the sign factor in 
A in (10 1 becomes a globally defined quantity (— 1)p a p b 



which agrees with the sign produced by commuting A 
and B when they share no common fermionic mode. 

The only assumption in the above presented fermionic 
rule is that one of the operators is of a well defined parity. 
Therefore all fermionic swap rules in RefsASHZEHI w ith 
more severe constraints (either both A and B are of well 
defined parity or even both are parity preserving) neces- 
sarily coincide with the fermionic rule presented here. 

Since the parity of Aij (and thus A'\ ■) is well defined 
by definition, one may use the swap rule to reverse the 
contraction order of A'] „• in ml such that 



A' 



1,1 



■A'i 



• • A' 



where operators A! { j are obtained from A t j by absorbing 
the local fermionic sign factor (— l)' +u arising from swap- 
ping two operators acting on a common virtual fermion, 



A i,3 - A lrudk\ i > a i, 3 l i+ 



-l,j a i,j+l c i,y 



(11) 



This way we are able to bring operators containing the 
same physical fermionic operator together and express 
the expectation value ^ of an arbitrary product opera- 
tor ([5]) in a form 

(*|O|*) = (-l)^^(0|^-«]...if[° 1 ^|0). (12) 

where operators K)??' 3 ^ are obtained by contracting over 
the physical mode as 



phys 



(13) 



In the language of tensor networks, this corresponds to 
obtaining a double-layer structure through contraction 
over the physical index in two single-layer structure of 
PEPS. Operators K^ 3 are of a well defined parity given 
by the corresponding operator Oi i.e. pij. Therefore, 
the contraction order may be chosen arbitrarily using the 
fermionic swap rule and anti commutation relations. In 
the following we will implicitly assume the dependence 
of K^fj''^ on the local operator Oij and use a compact 
notation Kij. 

Due to the canonical anti-commutation relations of 
fermionic operators, the tensor representation of K t j is 
not unique. Let us first choose a representation where 
the norm (^l^) is expressed as a tensor product 



) = tr^ 1 ' 11 ■ ■ • E [1 - n] ■ ■ ■ E [m - n] ) 



(14) 



where the multiplication order is given by the lattice 
bonds. Such form would enable us to contract the 
fermionic tensor network exactly. It is easy to show that 
this can be achieved by representing operators K^ j de- 
fined in ( 13 I in the following form 



Ki. 



i,3 " ^lrud a i,j + l^i+l,j a ij'Y; 



(15) 



where all (local) fermionic signs arising in the process are 
absorbed in tensor E^'fl. It should be noted, however, 
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that due to the fermionic signs the matrix r 'u'<i')«rud) 
is no longer positive semi-definite nor Hcrmitian as is 
the case in bosonic (spin) systems. Such an exact con- 
traction scheme is not limited to the calculation of the 
norm but can be used to contract exactly the expec- 
tation value of an arbitrary pro duct operator (J5|. In 
such a case, tensors E^'^ in ( 15 ) should be replaced by 

Efrld = E ff Li- 1 )- 12 J ' <3Pl+1 - 3 ' where PiJ is the P arit y 
of Oi.j in (j5f, defined globally. 

In order to draw the correspondence with PEPS algo- 
rithm, we will choose a different representation of Ki j 
where no signs are produced in contraction over a single 
row which will allow us to express the boundary row as a 
matrix product state. This is achieved by the following 
representation 



K 



i,j+lli+l,j^1i,j a i,j 



(16) 



where again all local signs are absorbed in tensor K_. 
For sake of concreteness, let us write the tensor elements 
Jf [i ' i] explicitly, 

KffJ = (-l) fKiL ^4l^ k ^A^ dk (17) 

with Offl = (0\c*' O u cl k \0) and sign function f K = VI + 
{I' + l)(r + u + d) + (r' +r)(u + d)+d(u' + u)+u'u + u' + 1' . 

The double layer structure given by pairs of fermionic 
operators a\ jOt'l j and similar, can also be interpreted as 

a structure given by higher-dimensional objects a\ j and 
similar. Evidently, the only property used in the formula- 
tion of fermionic network is the parity and all results also 
apply to higher-dimensional objects oq ^ where parity of 
cq is given by 



P(°k,j) = [p( a 'i,j) +p(a- J -)]mod2. 

The fact that the parity is the only relevant element 
in the anti-commutation relations of operators such as 
Ai j and i<Q j, suggests a natural way to generalize the 
tensor network to higher bond dimensions by replac- 
ing virtual fermionic operators ctij and jij in Q by 
higher dimensional objects ■ and 7. ., respectively. 

The rest of the method remains the same whereas all 
occupation numbers m appearing in fermionic sign fac- 
tors are replaced by the corresponding parities of m, i.e. 
p(m) = ( J^fc m k) mod 2. The only drawback in such gen- 
eralization scheme is that it confines the bond dimen- 
sion to the powers of two. An alternative generalization 
scheme is by combining fermionic and bosonic (spin) de- 
grees of freedom where the former would assure fermionic 
nature of description whereas the latter would enlarge the 
virtual space to capture more entangled physical systems. 
This would lead to bond dimensions that are even. 





FIG. 2. Merging two rows together ( 19 1 and replacing the 
double row by a MPS. 



B. Efficient contraction of fermionic tensor network 

Fermionic PEPS can be exactly contracted in a sign- 
free way using a suitable representation of operators Kij 
as shown in ( 14 ). However, exact contraction is only pos- 



sible with small bond dimension and small lattice sizes 
since the complexity scales exponentially with the lin- 
ear lattice dimension and we have to resort to an ap- 
proximate contraction scheme®. Let us quickly review 
the approximation scheme to calculate expectation val- 
ues as used in the PEPS algorithm. The first and the 
last rows are recognized as matrix product states £1 and 
£„, respectively, and all inner rows correspond to matrix 
product operators Ej for j G {2, 3, . . . , n — 1}. The ex- 
pectation value (^„|S n _iH rl _2 • • • S2 is calculated by 
approximating a product by a new matrix prod- 

uct state I £ 2 ) of some finite bond dimension^ and pro- 
ceeding iteratively until the expectation value is given by 

In the following we will show how the expectation value 
of an arbitrary product operator can be calculated effi- 
ciently in an approximate way, which is equivalent to 
the approach in occupation number representation^ or 
tensor network approach 12 ^. Taking the advantage of rep- 
resentation (16 1 of where no signs are produced in 
the horizontal contraction we express the first row as a 
matrix product state 



l2,n 



(18) 



[l,3,Oi,f] 
l.rfi.d 



The same applies 



with matrices (K [1j] -) Ll = K t 
to the last row. Inner rows, on the other hand, cannot be 
represented as matrix product operators in a form which 
would allow immediate contraction with matrix product 
states due to the fermionic signs produced by reordering 
vertical virtual fermionic operators. Nevertheless, using 
the fact that the parity of Kij is determined globally by 
the underlying operator Oij, one can change the con- 
traction order in contracting first two rows to 

■Ku = {-l) f K 2n K ln - ■ -K 21 K n (19) 



K2n- ■ -K 2 iKi n 
where / = Y^i= 



iPu^2j=iP'2j- Note that this step is 
trivial since there is no need for fermionic swap rules as 
no fermionic modes are crossed. Contracting products 
K 2 jK\j over the vertical mode (Fig. [2j and representing 



the result in form ( 16 1, we again obtain a matrix product 
description of form (18). This way, the fermionic nature 



is completely absorbed in the tensors and all the MPS 
formalism results apply. 
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C. Variational simulation of the ground state 

There are essentially two ways of simulating the ground 
state using tensor networks. The first possibility is the 
evolution of a PEPS state in imaginary time using the 
approximate Trotter decomposition of the evolution op- 
erator. The alternative way is to optimize tensors A^ 1 '^ 
site by site in a variational way such that the total energy 
E = (\&|i/|\f r )/(\f r |\I>) is minimal. While numerical stabil- 
ity often speaks in favor of the imaginary time evolution, 
the variational approach is faster and gives fairly good re- 
sults after a single optimization sweep over the lattice. In 
this paper, we shall only focus on the latter approach and 
show that all fermionic signs which appear in the com- 
putation are absorbed locally into tensors which makes 
the problem essentially sign-free for practical matters and 
thus well suited to conventional PEPS techniques. 

In the following we will show how to write the total 
energy as a function of a tensor A}*'^' in a sign-free way. 
Using the fermionic rule it is easy to show that the ex- 
pectation value (12 1 of an arbitrary product operator 
can be written as 



(*|o|*) = (o|IJ iJ .ng 1 I iJ .|o) 

where Ajj contains all tensor elements of A\ % ^ as 



(20) 



A, 



whereas flf^ contains tensors corresponding to all other 
lattice sites. Such expression may be easily obtained 
by replacing Kij in (12 1 according to (13 1 and anti- 



A [iJ] W r W d 

^ L l,r,u,d,k L 'i,j Lz i,j Lt i,j + l h,J h+l,j 



commuting Aij and Ajj to the far ends. Finally, all 
fermionic signs are absorbed in . Note that the sign 

in ( 12 ) is cancelled by commuting A! { j defined in ( fllj over 
all consequent sites and thus no longer appears m p0| . 
Using a convenient representation for fiffl, i.e. 



[O] 



n [0,i ' j] 

I'r'u'd'k'lrudk 
/f k' AX' t\ r' ,f w' ,/t d> 



(21) 



we are able to rewrite the expectation value as an ordi- 
nary scalar product 



(#|0|#> = A-nA 



(22) 



where vector elements of A are given by A^^} udk ^ and sim- 
ilarly for matrix elements of f2 given as £la' r 'u'd'k')(lrudk)- 
This way the expectation value of a fermionic operator 
is expressed in terms of a sign-free linear algebra ex- 
pression. Note however that the initial assumption that 

tensors A^ dk are of well defined parity, reduces the ef- 
fective subspace of the vector space to the even-even or 
odd-odd sector with respect to indices (I'r'u'd'k 1 ) and 
(Irudk). While operator il itself is always of even par- 
ity, i.e. p(l'r'u'd'k' Irudk) — 0, no such requirement is 



imposed separately to (I'r'u'd'k 1 ) and (Irudk). In princi- 
ple, both sub-sectors, even-even and odd-odd, should be 
obtained separately using the assumption for the parity 
of A or equivalently, tensor elements A^^ dk . However, 
since O is in total of even parity, no additional signs are 
produced in the odd-odd case where fl^ is represented 



in form ( 21 ) 



The total energy (*|if)|*)/(*|*) may be expressed 
in terms of effective operators as 



E = 



A ■ Keft A 
A-N cS A 



(23) 



where iV eff and EL e Q are obtained using the above de- 
scribed procedure for the identity operator and the 
Hamiltonian operator, respectively, where the latter is 
written as a superposition of product operators. Note 
that the computation of H_ c g is simplified for Hamilto- 
nians with local interactions where certain operators are 
grouped together in the (approximate) contraction pro- 
cess. 



The solution A which minimizes ( 23 ) is formally given 



by the lowest eigenvalue solution of a generalized eigen- 
value problem 



H cS A = XN oS A 



(24) 



where N c g is a semi-definite Hermitian matrix and iJ c g is 
Hermitian. Due to the parity constraints, the eigenvalue 
problem must be solved separately for both parity sub- 
sectors and the better solution should be retained. The 



generalized eigenvalue problem ( 24 ) is only well defined if 
iV e ff is nonsingular. In one-dimensional variational MPS 
with open boundary conditions one can always renor- 
malize the tensor network in a way that N c g is exactly 
equal to the identity which simplifies the computation 
and, more importantly, makes the method stable. In two 
dimensions, the way to make PEPS better conditioned re- 
mains an open question. In general, the spectrum of N c g 
might and does contain very small values or even zeros, in 
which case the standard algorithm would produce infinite 
or ill-disposed eigenvalues. The ill-conditioned general- 
ized eigenvalue problem must be solved in an approxi- 
mate fashion by isolating such invalid solutions either by 
projecting out the null-space of A[ Gff or using more so- 
phisticated algorithms such as Fix-Heiberger reduction^ 
where ill-conditioned modes of iV off are not completely 
neglected. 

We find that the most stable way is to project the sys- 
tem to the subspace spanned by well conditioned eigen- 
vectors of iV c ff with respect to a cutoff S and then using 
the Fix-Heiberger algorithm with the tolerance e > 105 
which eliminates all solutions unstable to the perturba- 
tion of e to the matrices -ff off and N_ g- In addition, 
when a good convergence is achieved, we optimize the 
total energy (23) in an iterative way using the conjugate- 
gradient method. Nevertheless, compromise between ef- 
ficiency and accuracy versus numerical stability must be 
made. 
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FIG. 3. Relative error of the ground state energy for the 
quadratic model [Eq. ??eq:Hl)] with 7 = 1 for lattice sizes 
10x10 and 4x4 and bond dimension D = 2 and D = 4. The 
truncation number is in all cases set to D = 64. 



III. RESULTS 



FIG. 4. Convergence in terms of the ground state energy 
relative error 5e (full lines) and the total particle number 
relative error <5jv (dotted lines) for ( 25 1 with 7 = 1 and A = 3. 
Three lattices sizes are considered as designated in the legend. 
Bond dimension is taken D — 4 with the truncation number 
D = 64. 



The finite size fermionic PEPS method is put to 
the test by simulating an exactly solvable bi-linear 
(quadratic) model on a square lattice. The model con- 
sists of three parts: hopping between nearest neighbor, 
pair creation/annihilation and chemical potential, de- 
scribed by the following Hamiltonian operator^ 



H 



[4( c * - 74)+h.-c- 



2^Ac£c„. 



(25) 



The pairing potential 7 > is used to destroy the to- 
tal particle number symmetry and A > is the chemical 
potential. The same model was also used irPSl where in- 
finite fermionic PEPS algorithm was presented. Unlike 
we assume open boundary conditions which is 
better suited for finite-size PEPS algorithm. The sys- 
tem is critical for A < 2 (gapless in the thermodynamic 
limit) and non-critical (gapped) elsewhere. We choose a 
line 7 = 1, A £ [1,3] in the parameter space and test 
the method with respect to the relative accuracy of the 
ground state energy as shown in Fig. [3J For the bond 
dimension we take either D = 2 or D — 4 which corre- 
sponds to one or two virtual fermions of each kind, re- 
spectively. The maximal bond dimension in the process 
of contracting the double-layer structure (see RefP for 
details) is designated by the truncation number D = 64 
which will be justified later. As expected, the method 
performs better in the gapped regime (A > 2) for both 
system sizes considered in Fig. [3j For the 10 x 10 lattice 
the spectral gap in the gapped regime at 7 = 1, A = 3 is 
of order of 3 • 10 3 1 £7 1 where E denotes the correspond- 
ing ground state energy and the total energy obtained 
from the simulation is below the energy of the first ex- 
cited state. In the gapless regime, e.g. for 7 = A = 1, no 
guarantee for the ground state is given since the spectral 
gap is of order of 10 _7 |i?o|- For the 4x4 lattice the spec- 



tral gap in the gapless regime is of order 10 3 1 £7o I which 
is a magnitude larger than the accuracy of the ground 
state energy. 

By increasing the bond dimension from D = 2 to 
D = 4, the relative accuracy is improved for an order 
of magnitude as seen in Fig. [3] for both lattice sizes 4x4 
and 10x10. Note however, that a fairly good precision is 
achieved already with the bond dimension D = 2. The al- 
gorithm would perform better for higher bond dimension 
if one could make PEPS well conditioned. Namely, with 



the increasing bond dimension the problem ( 23 1 becomes 



less conditioned and it is essential to use Fix-Heiberger 
procedure (and conjugate gradient method) to eliminate 
unstable solutions. If all nearly singular vectors were 
simply chopped away, the benefit of using higher bond 
dimension would be negligible. 

In Fig. [4] we show the convergence of the ground state 
energy and the total particle number as a function of the 
number K of single particle optimizations. The simula- 
tion is done first using the bond dimension D = 2 and 
switching to D = 4 when sufficiently good convergence 
rate (relative difference 10 -5 for the total energy between 
two consequent sweeps) is achieved. We consider three 
lattices sizes and observe that a fairly good approxima- 
tion to the ground state where the ground state energy 
is accurate to 1%, is achieved with less than two sweeps 
over the lattice. The initial state was in all cases taken 
random. After the initial sweep the convergence becomes 
slower but the relative error of both the ground state en- 
ergy and the total particle number typically decays as 
1/K. In variational methods such as PEPS the ground 
state energy is typically more accurate than other ob- 
servables such as the total particle number which is also 
confirmed in Fig. [4] We have however no explanation 
for the oscillations in accuracy for the total number of 
particles. 
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FIG. 5. Convergence for the ground state energy for ( 25 1 with 



■y=l, A = 2.5 on a 8x8 lattice. Legend entries designate 
D/D. Results for D — 4 are obtained starting from a good 
approximation for D = 2 (almost vertical lines on the main 
plot, magnified on the left in the linear-linear scale). 



Let us now check the validity of the results for various 
truncation numbers D used to truncate the large matrix 
products representing several consecutive rows. As pre- 
sented in Fig. [3J we used D = 64 which turned out to 
be sufficient to get good accuracy. In Fig. [5] we present 
the results for the quadratic model a 8x8 lattice where 
the same initial state was taken in all cases. We con- 
sider three different values of D for the bond dimension 
D = 2. Eventually, we start the simulation with the bond 
dimension D — 4 where the (almost converged) results 
from D = 2 were taken as the initial state, also magni- 
fied on the left side of Fig. [5J We observe that D — 32 is 
insufficient to achieve good accuracy of the ground state 
energy although it gives reasonable results with little ef- 
fort. There is virtually no difference between the cases 
D = 64 and D = 128 except the latter being computa- 
tionally much more demanding. As already mentioned in 
the previous section, the algorithm eventually produces 
unstable solutions where the effective norm operator A cff 
in p4| becomes more and more ill-conditioned. This is 
reflected in the oscillations seen in the magnification of 
Fig. [5] which are also a sign that the simulation should 
be stopped, unless the state is made better conditioned. 



IV. DISCUSSION 

The finite size fermionic PEPS method was tested 
on a trivial example of a quadratic integrable model 
where it was shown that fairly good results can be 
achieved for lattice sizes 10x10. The present formula- 
tion is however open to various improvements and modi- 
fications. The first improvement would be beneficial not 
just for fermionic PEPS but for all two dimensional PEPS 
Ansatze, namely a way to make PEPS better conditioned 
which is of crucial importance for employing higher bond 



FIG. 6. (Color online.) Convergence of the total energy Eq 
for the interacting model \2Q\ with V = 0.5 on a 4x4 lattice. 
Initial states are random with D — 4 (first two curves from 
the left), D = 2 (two curves on the right), and D — 2 switched 
to D = 4 (small deviation from the D — 2). The transition 
D = 2 to D = 4 is magnified in the plot. The inset shows the 
total particle number (JV) for the corresponding cases. The 
truncation number is set to D = 128 in all cases. 



dimensions. Another possibility would be to consider 
higher order symmetries such as the Zk symmetry for 
which the presented Zi symmetry algorithm presents a 
good starting point. 

The fermionic swapping rule allows arbitrary manip- 
ulations to the contraction order which enables various 
enhancements to the presented method. The first is a 
complementary way of optimizing tensors A} 1 '^ by imag- 
inary time evolution. The method can also be made more 
robust in convergence to the global minimum by adding 
stochastic updates to the tensor elements which would 
be beneficial especially with non-trivial models where the 
energy landscape is such that one easily gets stuck in a 
local minimum. Although such phenomenon was not ob- 
served in simulating the integrable model (25 1, it might 
occur for certain interacting models. Let us briefly con- 
sider an interacting model 



H 



■h.c] +V 



(26) 



where the total particle number (N) for N — JT ^ 
is a preserved quantity. The algorithm does not always 
converge to the global ground state but to the lowest- 
lying eigenstate in a particular total-particle number sub- 
sector, depending on the initial state. This issue may be 
addressed by simulating a modified model H' = H — [iN 
with the same eigenstates as ( 26 ) . Various total particle 



number sub-sectors are achieved by tuning the chemical 
potential. 

In Fig. [6] we present the total energy convergence for 
the interacting model (26) on a 4x4 lattice with the in- 
teraction strength V = 0.5. Different lines correspond to 
different random initial states and bond dimensions as 
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FIG. 7. (Color online.) Correlation functions (n,3,j«3,j'} (des- 
ignated as {rij-riji) in the graph, the second label is shown by 
dashed lines) for the interacting model ( |26[ ) with V = 0.5 on 
a 4x4 lattice. The data correspond to the left- most curve in 
the graph on Fig. [6] Exact values are plotted by thin gray 
lines. 



explained in the figure caption. We observe that in this 
case a quick convergence to the global minimum (note the 
exact energy levels designated by dotted lines) is achieved 
for all choices of random initial state although the rela- 
tive precision is not as good as in the integrable model. 
A higher bond dimension D = 4 gives better results after 
fewer number of iteration steps but effectively consumes 
more computational time. Similar accuracy is obtained 
when an approximate ground state is obtained by a small 
bond dimension D = 2 and later switched to D — 4. The 
simulation however quickly stops due to achieved relative 
accuracy between subsequent sweeps. The fluctuations 
in the energy in D = 4 are explained by the transitions 
between rows when an error is made in truncating large 
matrix product states (note that the energy is calculated 
in an approximate way). The total particle number is in 
all cases in agreement with the exact value in the ground 
state, (N) = 6 for this choice of parameters. It must be 
noted that not every initial state converges to the ground 
state but might as well converge to a local minimum. No 
such case was however observed for V — 0.5 on a 4x4 
lattice. It might be beneficial to tune the chemical po- 
tential to influence the number of particles in the system 
or start with a good initial state. 

Apart from observables consisting of local contribu- 



tions such as the energy or the total particle number, 
we can also investigate nonlocal quantities such as cor- 
relation functions, e.g. the density-density correlations 
(riijnii ji) . In Fig.[7]we show the density-density correla- 
tions for a fixed row (i = 3) for the interacting model ( 26 ) 
on a 4x4 lattice. The data correspond to the first curve 
from the left in Fig. [6] For computational simplicity we 
only calculate the correlations after each complete sweep, 
i.e. every (2mn — 2)th step. It can be seen that in all 
cases the results practically coincide with the exact re- 
sults designated by thin gray horizontal lines with the 
absolute error of order of 10 -3 as shown in the inset. 



V. CONCLUSION 

We have presented a finite size fermionic PEPS method 
to simulate ground states of two dimensional fermionic 
systems-^ completely in terms of fermionic operators. 
Using a fermionic swap rule to reverse the contraction or- 
der of two superpositions of products of fermionic canon- 
ical operators we have shown how a fermionic tensor net- 
work can be contracted exactly without introducing any 
additional sign bond but instead absorbing all signs lo- 
cally. Due to the parity constraints in fermionic systems 
we have presented a way, equivalent td 17 * 2 ° l, to calculate 
the expectation values for arbitrary operators efficiently 
in an approximate fashion. Finally, we have implemented 
the variational PEPS algorithm on a fermionic lattice and 
tested it on an integrable bi-linear fermionic model. We 
have found that the ground states of such a model can be 
simulated efficiently with relatively high accuracy in the 
ground state energy and the total number of particles. 
We have also discussed the performance of the method 
in the case of an interacting fermionic system where the 
method converges to the global minimum, albeit with 
less accurate precision. Besides local observables such as 
the energy and the total number of particles, the method 
correctly describes also the non-local two-point correla- 
tions. 
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